{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hello Radiomics example: using the feature extractor to calculate features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example shows how to use the radiomics package and the feature extractor.\n",
    "The feature extractor handles preprocessing, and then calls the needed featureclasses to calculate the features.\n",
    "It is also possible to directly instantiate the feature classes. However, this is not recommended for use outside debugging or development. For more information, see `helloFeatureClass`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import annotations\n",
    "\n",
    "import logging\n",
    "import os\n",
    "\n",
    "import radiomics\n",
    "from radiomics import featureextractor, getFeatureClasses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up logging"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regulate verbosity of PyRadiomics (outputs to stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Regulate verbosity with radiomics.setVerbosity\n",
    "# radiomics.setVerbosity(logging.INFO)  # Use logging.DEBUG for maximum output, default verbosity level = WARNING"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up logging to a log file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Get the PyRadiomics logger (default log-level = INFO)\n",
    "logger = radiomics.logger\n",
    "logger.setLevel(logging.DEBUG)  # set level to DEBUG to include debug log messages in log file\n",
    "\n",
    "# Write out all log entries to a file\n",
    "handler = logging.FileHandler(filename='testLog.txt', mode='w')\n",
    "formatter = logging.Formatter('%(levelname)s:%(name)s: %(message)s')\n",
    "handler.setFormatter(formatter)\n",
    "logger.addHandler(handler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting the testcase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test cases can be downloaded to temporary files. This is handled by the `radiomics.getTestCase()` function, which checks if the requested test case is available and if not, downloads it. It returns a tuple with the location of the image and mask of the requested test case, or (None, None) if it fails.\n",
    "\n",
    "Alternatively, if the data is available somewhere locally, this directory can be passed as a second argument to `radiomics.getTestCase()`. If that directory does not exist or does not contain the testcase, functionality reverts to default and tries to download the test data.\n",
    "\n",
    "If getting the test case fails, PyRadiomics will log an error explaining the cause."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "featureClasses = getFeatureClasses()\n",
    "imageName, maskName = radiomics.getTestCase('brain1')\n",
    "\n",
    "if imageName is None or maskName is None:  # Something went wrong, in this case PyRadiomics will also log an error\n",
    "    msg = 'Error getting testcase!'\n",
    "    raise Exception(msg)  # Raise exception to prevent cells below from running in case of \"run all\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initializing the feature extractor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extraction Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use a parameter file, this customizes the extraction settings and also specifies the input image types to use and which features should be extracted\n",
    "params = os.path.join('..', 'examples', 'exampleSettings', 'Params.yaml')\n",
    "\n",
    "extractor = featureextractor.RadiomicsFeatureExtractor(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Alternative: use hardcoded settings (separate for settings, input image types and enabled features)\n",
    "settings = {}\n",
    "settings['binWidth'] = 25\n",
    "settings['resampledPixelSpacing'] = None\n",
    "# settings['resampledPixelSpacing'] = [3, 3, 3]  # This is an example for defining resampling (voxels with size 3x3x3mm)\n",
    "settings['interpolator'] = 'sitkBSpline'\n",
    "settings['verbose'] = True\n",
    "\n",
    "extractor = featureextractor.RadiomicsFeatureExtractor(**settings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Input images: applying filters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enabled input images:\n",
      "\tOriginal\n"
     ]
    }
   ],
   "source": [
    "# By default, only 'Original' (no filter applied) is enabled. Optionally enable some image types:\n",
    "\n",
    "# extractor.enableImageTypeByName('Wavelet')\n",
    "# extractor.enableImageTypeByName('LoG', customArgs={'sigma':[3.0]})\n",
    "# extractor.enableImageTypeByName('Square')\n",
    "# extractor.enableImageTypeByName('SquareRoot')\n",
    "# extractor.enableImageTypeByName('Exponential')\n",
    "# extractor.enableImageTypeByName('Logarithm')\n",
    "\n",
    "# Alternative; set filters in one operation \n",
    "# This updates current enabled image types, i.e. overwrites custom settings specified per filter. \n",
    "# However, image types already enabled, but not passed in this call, are not disabled or altered.\n",
    "\n",
    "# extractor.enableImageTypes(Wavelet={}, LoG={'sigma':[3.0]})\n",
    "\n",
    "print('Enabled input images:')\n",
    "for imageType in extractor.enabledImagetypes:\n",
    "    print('\\t' + imageType)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Feature classes: setting which feature(classes) need to be calculated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Disable all classes\n",
    "extractor.disableAllFeatures()\n",
    "\n",
    "# Enable all features in firstorder\n",
    "extractor.enableFeatureClassByName('firstorder')\n",
    "\n",
    "# Alternative; only enable 'Mean' and 'Skewness' features in firstorder\n",
    "# extractor.enableFeaturesByName(firstorder=['Mean', 'Skewness'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting the docstrings of the active features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Active features:\n",
      "InterquartileRange\n",
      "\n",
      "    **10. Interquartile Range**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{interquartile range} = \\textbf{P}_{75} - \\textbf{P}_{25}\n",
      "\n",
      "    Here :math:`\\textbf{P}_{25}` and :math:`\\textbf{P}_{75}` are the 25\\ :sup:`th` and 75\\ :sup:`th` percentile of the\n",
      "    image array, respectively.\n",
      "    \n",
      "Skewness\n",
      "\n",
      "    **16. Skewness**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{skewness} = \\displaystyle\\frac{\\mu_3}{\\sigma^3} =\n",
      "      \\frac{\\frac{1}{N_p}\\sum^{N_p}_{i=1}{(\\textbf{X}(i)-\\bar{X})^3}}\n",
      "      {\\left(\\sqrt{\\frac{1}{N_p}\\sum^{N_p}_{i=1}{(\\textbf{X}(i)-\\bar{X})^2}}\\right)^3}\n",
      "\n",
      "    Where :math:`\\mu_3` is the 3\\ :sup:`rd` central moment.\n",
      "\n",
      "    Skewness measures the asymmetry of the distribution of values about the Mean value. Depending on where the tail is\n",
      "    elongated and the mass of the distribution is concentrated, this value can be positive or negative.\n",
      "\n",
      "    Related links:\n",
      "\n",
      "    https://en.wikipedia.org/wiki/Skewness\n",
      "\n",
      "    .. note::\n",
      "      In case of a flat region, the standard deviation and 4\\ :sup:`rd` central moment will be both 0. In this case, a\n",
      "      value of 0 is returned.\n",
      "    \n",
      "TotalEnergy\n",
      "\n",
      "    **2. Total Energy**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{total energy} = V_{voxel}\\displaystyle\\sum^{N_p}_{i=1}{(\\textbf{X}(i) + c)^2}\n",
      "\n",
      "    Here, :math:`c` is optional value, defined by ``voxelArrayShift``, which shifts the intensities to prevent negative\n",
      "    values in :math:`\\textbf{X}`. This ensures that voxels with the lowest gray values contribute the least to Energy,\n",
      "    instead of voxels with gray level intensity closest to 0.\n",
      "\n",
      "    Total Energy is the value of Energy feature scaled by the volume of the voxel in cubic mm.\n",
      "\n",
      "    .. note::\n",
      "      This feature is volume-confounded, a larger value of :math:`c` increases the effect of volume-confounding.\n",
      "\n",
      "    .. note::\n",
      "      Not present in IBSI feature definitions\n",
      "    \n",
      "Uniformity\n",
      "\n",
      "    **19. Uniformity**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{uniformity} = \\displaystyle\\sum^{N_g}_{i=1}{p(i)^2}\n",
      "\n",
      "    Uniformity is a measure of the sum of the squares of each intensity value. This is a measure of the homogeneity of\n",
      "    the image array, where a greater uniformity implies a greater homogeneity or a smaller range of discrete intensity\n",
      "    values.\n",
      "\n",
      "    .. note::\n",
      "      Defined by IBSI as Intensity Histogram Uniformity.\n",
      "    \n",
      "Median\n",
      "\n",
      "    **9. Median**\n",
      "\n",
      "    The median gray level intensity within the ROI.\n",
      "    \n",
      "Energy\n",
      "\n",
      "    **1. Energy**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{energy} = \\displaystyle\\sum^{N_p}_{i=1}{(\\textbf{X}(i) + c)^2}\n",
      "\n",
      "    Here, :math:`c` is optional value, defined by ``voxelArrayShift``, which shifts the intensities to prevent negative\n",
      "    values in :math:`\\textbf{X}`. This ensures that voxels with the lowest gray values contribute the least to Energy,\n",
      "    instead of voxels with gray level intensity closest to 0.\n",
      "\n",
      "    Energy is a measure of the magnitude of voxel values in an image. A larger values implies a greater sum of the\n",
      "    squares of these values.\n",
      "\n",
      "    .. note::\n",
      "      This feature is volume-confounded, a larger value of :math:`c` increases the effect of volume-confounding.\n",
      "    \n",
      "RobustMeanAbsoluteDeviation\n",
      "\n",
      "    **13. Robust Mean Absolute Deviation (rMAD)**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{rMAD} = \\frac{1}{N_{10-90}}\\displaystyle\\sum^{N_{10-90}}_{i=1}\n",
      "      {|\\textbf{X}_{10-90}(i)-\\bar{X}_{10-90}|}\n",
      "\n",
      "    Robust Mean Absolute Deviation is the mean distance of all intensity values\n",
      "    from the Mean Value calculated on the subset of image array with gray levels in between, or equal\n",
      "    to the 10\\ :sup:`th` and 90\\ :sup:`th` percentile.\n",
      "    \n",
      "MeanAbsoluteDeviation\n",
      "\n",
      "    **12. Mean Absolute Deviation (MAD)**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{MAD} = \\frac{1}{N_p}\\displaystyle\\sum^{N_p}_{i=1}{|\\textbf{X}(i)-\\bar{X}|}\n",
      "\n",
      "    Mean Absolute Deviation is the mean distance of all intensity values from the Mean Value of the image array.\n",
      "    \n",
      "Maximum\n",
      "\n",
      "    **7. Maximum**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{maximum} = \\max(\\textbf{X})\n",
      "\n",
      "    The maximum gray level intensity within the ROI.\n",
      "    \n",
      "RootMeanSquared\n",
      "\n",
      "    **14. Root Mean Squared (RMS)**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{RMS} = \\sqrt{\\frac{1}{N_p}\\sum^{N_p}_{i=1}{(\\textbf{X}(i) + c)^2}}\n",
      "\n",
      "    Here, :math:`c` is optional value, defined by ``voxelArrayShift``, which shifts the intensities to prevent negative\n",
      "    values in :math:`\\textbf{X}`. This ensures that voxels with the lowest gray values contribute the least to RMS,\n",
      "    instead of voxels with gray level intensity closest to 0.\n",
      "\n",
      "    RMS is the square-root of the mean of all the squared intensity values. It is another measure of the magnitude of\n",
      "    the image values. This feature is volume-confounded, a larger value of :math:`c` increases the effect of\n",
      "    volume-confounding.\n",
      "    \n",
      "90Percentile\n",
      "\n",
      "    **6. 90th percentile**\n",
      "\n",
      "    The 90\\ :sup:`th` percentile of :math:`\\textbf{X}`\n",
      "    \n",
      "Minimum\n",
      "\n",
      "    **4. Minimum**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{minimum} = \\min(\\textbf{X})\n",
      "    \n",
      "Entropy\n",
      "\n",
      "    **3. Entropy**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{entropy} = -\\displaystyle\\sum^{N_g}_{i=1}{p(i)\\log_2\\big(p(i)+\\epsilon\\big)}\n",
      "\n",
      "    Here, :math:`\\epsilon` is an arbitrarily small positive number (:math:`\\approx 2.2\\times10^{-16}`).\n",
      "\n",
      "    Entropy specifies the uncertainty/randomness in the image values. It measures the average amount of information\n",
      "    required to encode the image values.\n",
      "\n",
      "    .. note::\n",
      "      Defined by IBSI as Intensity Histogram Entropy.\n",
      "    \n",
      "Range\n",
      "\n",
      "    **11. Range**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{range} = \\max(\\textbf{X}) - \\min(\\textbf{X})\n",
      "\n",
      "    The range of gray values in the ROI.\n",
      "    \n",
      "Variance\n",
      "\n",
      "    **18. Variance**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{variance} = \\frac{1}{N_p}\\displaystyle\\sum^{N_p}_{i=1}{(\\textbf{X}(i)-\\bar{X})^2}\n",
      "\n",
      "    Variance is the the mean of the squared distances of each intensity value from the Mean value. This is a measure of\n",
      "    the spread of the distribution about the mean. By definition, :math:`\\textit{variance} = \\sigma^2`\n",
      "    \n",
      "10Percentile\n",
      "\n",
      "    **5. 10th percentile**\n",
      "\n",
      "    The 10\\ :sup:`th` percentile of :math:`\\textbf{X}`\n",
      "    \n",
      "Kurtosis\n",
      "\n",
      "    **17. Kurtosis**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{kurtosis} = \\displaystyle\\frac{\\mu_4}{\\sigma^4} =\n",
      "      \\frac{\\frac{1}{N_p}\\sum^{N_p}_{i=1}{(\\textbf{X}(i)-\\bar{X})^4}}\n",
      "      {\\left(\\frac{1}{N_p}\\sum^{N_p}_{i=1}{(\\textbf{X}(i)-\\bar{X}})^2\\right)^2}\n",
      "\n",
      "    Where :math:`\\mu_4` is the 4\\ :sup:`th` central moment.\n",
      "\n",
      "    Kurtosis is a measure of the 'peakedness' of the distribution of values in the image ROI. A higher kurtosis implies\n",
      "    that the mass of the distribution is concentrated towards the tail(s) rather than towards the mean. A lower kurtosis\n",
      "    implies the reverse: that the mass of the distribution is concentrated towards a spike near the Mean value.\n",
      "\n",
      "    Related links:\n",
      "\n",
      "    https://en.wikipedia.org/wiki/Kurtosis\n",
      "\n",
      "    .. note::\n",
      "      In case of a flat region, the standard deviation and 4\\ :sup:`rd` central moment will be both 0. In this case, a\n",
      "      value of 0 is returned.\n",
      "\n",
      "    .. note::\n",
      "      The IBSI feature definition implements excess kurtosis, where kurtosis is corrected by -3, yielding 0 for normal\n",
      "      distributions. The PyRadiomics kurtosis is not corrected, yielding a value 3 higher than the IBSI kurtosis.\n",
      "    \n",
      "Mean\n",
      "\n",
      "    **8. Mean**\n",
      "\n",
      "    .. math::\n",
      "      \\textit{mean} = \\frac{1}{N_p}\\displaystyle\\sum^{N_p}_{i=1}{\\textbf{X}(i)}\n",
      "\n",
      "    The average gray level intensity within the ROI.\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "print('Active features:')\n",
    "for cls, orig_features in extractor.enabledFeatures.items():\n",
    "    features = orig_features\n",
    "    if len(features) == 0:\n",
    "        features = [f for f, deprecated in featureClasses[cls].getFeatureNames().items() if not deprecated]\n",
    "    for f in features:\n",
    "        print(f)\n",
    "        print(getattr(featureClasses[cls], f'get{f}FeatureValue').__doc__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the values of the active features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating features\n"
     ]
    }
   ],
   "source": [
    "print('Calculating features')\n",
    "featureVector = extractor.execute(imageName, maskName)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed diagnostics_Versions_PyRadiomics: 2.1.2.post68.dev0+g590d9fe\n",
      "Computed diagnostics_Versions_Numpy: 1.15.1\n",
      "Computed diagnostics_Versions_SimpleITK: 0.9.1\n",
      "Computed diagnostics_Versions_PyWavelet: 0.5.2\n",
      "Computed diagnostics_Versions_Python: 2.7.11\n",
      "Computed diagnostics_Configuration_Settings: {'distances': [1], 'verbose': True, 'additionalInfo': True, 'force2D': False, 'interpolator': 'sitkBSpline', 'resampledPixelSpacing': None, 'label': 1, 'normalizeScale': 1, 'normalize': False, 'force2Ddimension': 0, 'removeOutliers': None, 'minimumROISize': None, 'binWidth': 25, 'minimumROIDimensions': 2, 'preCrop': False, 'resegmentRange': None, 'padDistance': 5}\n",
      "Computed diagnostics_Configuration_EnabledImageTypes: {'Original': {}}\n",
      "Computed diagnostics_Image-original_Hash: 5c9ce3ca174f0f8324aa4d277e0fef82dc5ac566\n",
      "Computed diagnostics_Image-original_Dimensionality: 3D\n",
      "Computed diagnostics_Image-original_Spacing: (0.7812499999999999, 0.7812499999999999, 6.499999999999998)\n",
      "Computed diagnostics_Image-original_Size: (256, 256, 25)\n",
      "Computed diagnostics_Image-original_Mean: 385.6564080810547\n",
      "Computed diagnostics_Image-original_Minimum: 0.0\n",
      "Computed diagnostics_Image-original_Maximum: 3057.0\n",
      "Computed diagnostics_Mask-original_Hash: 9dc2c3137b31fd872997d92c9a92d5178126d9d3\n",
      "Computed diagnostics_Mask-original_Spacing: (0.7812499999999999, 0.7812499999999999, 6.499999999999998)\n",
      "Computed diagnostics_Mask-original_Size: (256, 256, 25)\n",
      "Computed diagnostics_Mask-original_BoundingBox: (162, 84, 11, 47, 70, 7)\n",
      "Computed diagnostics_Mask-original_VoxelNum: 4137\n",
      "Computed diagnostics_Mask-original_VolumeNum: 2\n",
      "Computed diagnostics_Mask-original_CenterOfMassIndex: (186.98549673676578, 106.3562968334542, 14.38917089678511)\n",
      "Computed diagnostics_Mask-original_CenterOfMass: (46.47304432559825, 16.518518098863908, 15.529610829103234)\n",
      "Computed original_firstorder_InterquartileRange: 253.0\n",
      "Computed original_firstorder_Skewness: 0.27565085908587594\n",
      "Computed original_firstorder_Uniformity: 0.045156963555862184\n",
      "Computed original_firstorder_Median: 812.0\n",
      "Computed original_firstorder_Energy: 2918821481.0\n",
      "Computed original_firstorder_RobustMeanAbsoluteDeviation: 103.00138343026683\n",
      "Computed original_firstorder_MeanAbsoluteDeviation: 133.44726195252767\n",
      "Computed original_firstorder_TotalEnergy: 11579797135.314934\n",
      "Computed original_firstorder_Maximum: 1266.0\n",
      "Computed original_firstorder_RootMeanSquared: 839.9646448180755\n",
      "Computed original_firstorder_90Percentile: 1044.4\n",
      "Computed original_firstorder_Minimum: 468.0\n",
      "Computed original_firstorder_Entropy: 4.601935553903786\n",
      "Computed original_firstorder_Range: 798.0\n",
      "Computed original_firstorder_Variance: 24527.07920837261\n",
      "Computed original_firstorder_10Percentile: 632.0\n",
      "Computed original_firstorder_Kurtosis: 2.1807729393860265\n",
      "Computed original_firstorder_Mean: 825.2354363065023\n"
     ]
    }
   ],
   "source": [
    "# Show output\n",
    "for featureName in featureVector:\n",
    "    print(f'Computed {featureName}: {featureVector[featureName]}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
